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Abstract 

EVITA, standing for Evolutionary Inventory and Transportation Al- 
gorithm, is a two-level methodology designed to address the Inventory 
and Transportation Problem (ITP) in retail chains. The top level uses 
an evolutionary algorithm to obtain delivery patterns for each shop on 
a weekly basis so as to minimise the inventory costs, while the bottom 
level solves the Vehicle Routing Problem (VRP) for every day in order 
to obtain the minimum transport costs associated to a particular set of 
patterns. 

The aim of this paper is to investigate whether a multiobjective ap- 
proach to this problem can yield any advantage over the previously used 
single objective approach. The analysis performed allows us to conclude 
that this is not the case and that the single objective approach is in gene- 
ral preferable for the ITP in the case studied. A further conclusion is that 
it is useful to employ a classical algorithm such as Clarke & Wright's as 
the seed for other metaheuristics like local search or tabu search in order 
to provide good results for the Vehicle Routing Problem. 



1 Introduction 

Given a retail chain and a central depot that supplies it, both belonging to the 
same company, we define the Inventory and Transportation Problem as that 
whose objective is to minimise the costs of both inventory and transportation, 
subject to a number of constraints imposed at the shop level. 
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In previous work ( Esparcia- Alcazar et al. . 2006al lb. 2007af bl. 20091 ) we em- 



ployed a single objective approach to this problem, which aimed at minimising 
the total weekly cost calculated as the sum of the inventory and transportation 
costs. The purpose of the current work is to evaluate the convenience of adopt- 
ing a multiobjective approach to the problem. Although it is true that the aim 
of the ITP is to minimise the global cost and this is calculated by simply adding 
the costs of inventory and transport (both measured in currency units), it is no 
less true that the minimisation of both costs are contradictory objectives since, 
as was stated above, reducing the cost of transport implies increasing the cost 
of inventory and vice versa. For this reason it is possible that a multiobjective 
approach will obtain better results in this problem, and this is what we set out 
to investigate. 

With this aim in mind we have carried out an extensive series of experiments 
using the eight instances used in ( Esparcia- Alcazar et al. . 20091 ). plus two new 
oncflj. The set of restrictions, consisting of the characteristics of the vehicles 
employed and the working hours of the drivers (see Table [7]), plus the parameter 
configuration of the evolutionary algorithm are kept as in that work. 



In the multiobjective approach we have used the NSGA-II algorithm (jPeb et al 

2000h . which is described in more detail in Appendix I. 



A further objective is related to the fact that the choice of the algorithm that 
yields the transportation routes (the VR P solver) can play a significant part in 



the performance of the whole algorithm (jEsparcia- Alcazar et al.l . l2006af ). Here 



we will employ improved versions of three different algorithms that have been 
employed in the literature for solving the VRP and will compare them on the se- 
lected problem instances, with the aim to determine whether there is a relation 
between the instance characteristics and the performance of the VRP solver. 
The VRP algo rithms studied are improved variants of th e three that obtained 



best results in (jEsparcia- Alcazar et al.l.l2009): tab u search (jCordeau et al.l . ll997f) 



ant colony optimisation ( Dong and XianJ 



20061) and a classical VR P solving 



technique, Clarke and Wright's algorithm (jClarke and Wrightl . 119641 ). 

The rest of the paper is structured as follows. Section [2] provides background 
on the problem; it contains a summary of the state of the art in this and related 
problems and a detailed description of the ITP. The top and lower levels of the 
algorithm are described in Sections [3] and [H the latter containing the details of 
the three VRP solvers employed. The experimental setup is given in Section 
[5l with results and analysis thereof contained in Section \6\ Finally, Section [7] 
presents the conclusions and outlines future areas of research. 



2 Problem background 

The ITP deals with the management of two different aspects of retail chain 
logistics, inventory costs and transportation costs, which have both received 
lots of attention in the logistics literature. Here we will describe how they have 

^Nine instances are t aken from |http: //branchandcut . org/VRP/data/| and the remaining 
one from |http : //www. f ermml-hagen.de/WIMF/touren/lrLhalte/problnst .htm| 
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been addressed in the past and what aspects are particularly relevant to our 
case. 



2.1 The Vehicle Routing Problem 

The Vehicle Routing Problem (VRP) consists on finding an optimal set of del i 
very routes from a depot to a set of customers to serve ( Toth and Vigol 2001 ) 



The routes must start and finish in the depot and each customer must be served 
by one and only one vehicle^, which means a customer cannot be contained in 
more than one route. Different versions of the problem have slightly different 
objectives or ways to define optimality: it can refer to finding the minimum 
cost, employing minimum time, minimum number of delivery vehicles or com- 
binations of these and other factors. 

Amongst these variants of the problem the most popular is the capacitated 
VRP, or CVRP, which refers to the fact that each delivery vehicle has a limited 
capacit}0. Also popular is the VRP with time windows, or VRPTW, in which 
each customer must be served during a specified time interval or window. 

Another variant of the VRP relevant to our problem is the periodic vehicle 
routing problem (PVRP), which appears when customers have established a 
predetermined delivery frequency and a combination of admissible delivery days 
within the planning horizon. The objective is to minimise the total duration 
of the routes, while the restrictions usually involve a limited capacity of the 
delivery vehicle s and a max imum duration of eac h itinerary. See for instance, 
( Cordeau et al. . 1997 ) and ( Toth and Vigol . 2001) for a general description of 



the periodic VRP. 

In this paper we are concerned with the CVRP, simply referred to as VRP. 
The only limitation will be in the capacity of the vehicles used, and not in their 
number. We will also consider that the fleet is homogeneous, i.e. only one type 
of vehicle is used, with unique values of average velocity and capacity. 

For simplicity reasons we will consider that the customers (shops in the retail 
chain in our case) have no time windows, i.e. the deliveries can take place at any 
time. However, the hours a driver can work are limited by regulations and this 
has to be taken into account in the time needed per delivery, as is the unloading 
time. Also, it should be possible to serve all customers; this means that routes 
from the depot to any customer and back must take less than the working hours 
of the driver. 

Finally, we will consider that the transport cost only includes the cost per 
kilometer; there is no cost attached to either the time of use of the vehicle or to 
the number of units delivered, nor there is a fixed cost per vehicle. 



^No single customer can consume more than the capacity of any one vehicle. 

VRP in which vehicles had infinite capacity could be subsumed under the general 
Traveling Salesman Problem 
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2.2 Inventory and transportation management 

One step beyond the VRP is the pro blem of optimising simul t aneous ly the costs 
of inventory and transportation. In ( Esparcia- Alcazar et all [2007al ) the reader 
can find a review of different approaches employed in the literature. 

Amongst these, one of the most relevant to our case is the inventory rou- 
ting problem (IRP), which arises when a vendor del ivers a single product and 
implements a Vendor Inventory Management (VIM) ( (^etinkava and Le3 . 2000f) 
policy with its clients, so that the vendor decides the delivery (time and quan- 
tity) in order to prevent the clients from run ning out of stock while min i misin g 
transportation and inventory holding costs (jCampbell and Savelsberghl . 120041 ) . 
However, retail chains cannot be addressed in this way, as thousands of items 
are involved. 

The work presented here focuses on the Inventory and Transportation Pro - 



200d) 



blem (ITP), which was first described in (jCardos and Garcfa-Sabater 
with the aim of addressing the case of retail chains. Thus, the ITP can be 
viewed as a generalisation of the IRP to the multiproduct case. Additionally, it 
can also be viewed as a variant of the PVRP described previously that includes 
inventory cost^ and a set of delivery frequencies instead of a unique delivery 
frequency for each shop. 

The main feature that differentiates the ITP f rom other similar s upply 

chain management ones addressed in the literature dCetinkava and Lee ^ 2000 ; 

dos Santos Coelho and Loped. 2006 : Federgruen and Zipkin . 1984: Sin dhuchao et al 
200,4 IViswanathan and Mathuii . ll997[ ) is that we have to decide on the frequency 
of delivery to each shop, which determines the size of the deliveries. The inven- 
tory costs can then be calculated accordingly, assuming a commonly-employed 
periodic review stock policy for the retail chain shops. Besides, for a given de- 
livery frequency, expressed in terms of number of days a week, there can also 
be a number of delivery patterns, i.e. the specific days of the week in which 
the shop is served. Once these are established, the transportation costs can be 
calculated by solving the VRP for each day of the week. 

Because a pattern assumes a given frequency, the problem is limited to 
obtaining the optimal patterns (one per shop) and set of routes (one set per day) . 
The optimum is defined as a combination of patterns and routes that minimises 
the total cost, which is calculated as the sum of the individual inventory costs 
per shop (inventory cost) plus the sum of the transportation costs for all days 
of the week (transport cost). These two objectives are in general contradictory: 
the higher the frequency of delivery the lower the inventory cost, but conversely, 
a higher frequency involves higher transportation costs. 

The operational constraints a t the shop level are imposed by the business 
logic and can be listed as follows ( Esparcia- Alcazar et al. . 20091 ): 



1. A periodic review stock policy is applied for the shop items. This means 
that the decision of whether to deliver to a particular shop is taken cen- 
trally and not at the shop. As a consequence, stockout is allowed. 



I.e. the PVRP can be seen as an ITP in which the inventory costs are zero. 
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2. Shops have a hmited stock capacity. 

3. The retail chain tries to fulfil backorders in as few days as possible, so 
there is a lower bound for delivery frequency, which depends on the target 
client service level. 

4. The expected stock reduction between replenishments (deliveries) cannot 
be too high in order to avoid two problems: (a) the unappealing empty- 
shelves aspect of the shop just before the replenishment; and (b) replen- 
ishment orders too large to be placed on the shelves by the shop personnel 
in a short time compatible with their primary selling activity. 

5. Conversely, the expected stock reduction must be high enough to perform 
an efficient allocation of the replenishment order. 

6. As a consequence of points [2] to [5] above, not all frequencies are admissible 
for all shops; in general, very high and very low frequencies are undesirable. 
For instance, frequency 1 (one delivery per week) is not applicable to most 
shops. 

7. Sales are not uniformly distributed over the time horizon (week), ten- 
ding to increase over the weekend. Hence, in order to match deliveries 
to sales, only a given number of delivery patterns are allowed for every 
feasible frequency. For instance, a frequency-2 pattern such as (Mon, Fri) 
is admissible, while another of the same frequency such as (Mon, Tues) is 
not. 

8. Although we are dealing with thousands of items, the load is containerised; 
hence, the size of the deliveries is expressed as an integer, representing the 
number of roll-containers. 

Thus, to summarise, our task involves finding: 

• The optimal set of patterns, Vopt, with which all shops can be served. A 
pattern p represents a set of days in which a shop is served which implies 
a delivery frequency (expressed as number of days a week) for the shop 

• The optimal routes for each day of the working week, by solving the VRP 
for the shops allocated to that day by the corresponding pattern 

2.3 Objective functions in EVITA 

EVITA operates in two levels: the lower one deals exclusively with the trans- 
portation costs per day and the top one incorporates these and the inventory 
costs into the total costs. 

In the single objective case, the optimum is defined as the solution that 
minimises the total cost, given by the function 

/ = TotalCost = InventoryC ost + TransportCost (1) 
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In the multiple objective approach, we will deal with two separate cost func- 
tions: 



fi — InventoryCost (2) 
ft = TransportCost (3) 

However, we will still use the total cost defined for the single objective case 
to compare the results between multi and monoobjective solutions. The reason 
is that the company is interested in spending less, irrespective of where the 
reduction comes from, and, at the end of the day, both costs come in euros. 

Inventory costs are computed from the patterns for each shop by taking into 
account the associated delivery frequency and looking up the inventory cost per 
shop in the corresponding table. An example of the latter is given in Table [TJ 



Table 1: Inventory cost (in Euro) and size of the deliveries per shop (expressed 
in roll containers) depending on the delivery frequency (in days). Missing data 
corresponds to frequencies that are not admissible for each shop. 





Inventory cost 


Delivery size 






(€) 


(roll containers) 




Shop # 


Frequency (days) 


Frequency (days) 






1 2 3 4 5 


12 3 4 


5 


1 


- - - 336 325 


- - - 2 


2 


2 


- - - 335 325 


- - - 2 


2 


N 


- 311 293 286 284 


- 3 2 2 


1 



For instance, let us assume that shop N was assigned a pattern of frequency 
4; we would look up in the table the inventory cost for the shop at that frequency, 
which is 286€. Proceeding in the same way with all shops and adding up the 
results we would obtain the total inventory cost. 

Transport costs are obtained by solving the VRP with one of the algorithms 
under study. The demands (size of deliveries) of each shop would also be taken 
from Table [H In the example above, for shop N at frequency 4 the delivery size 
is 2 roll containers. 

Problem data is freely available from our group website: http : / / casnew . it i . es^ 



^Its use is subject to the condition that this or other papers on the same subject by the 
authors are mentioned 
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3 The top level: Evolutionary Algorithm 



The top level in EVITA is an evolutionary algorithm in which a population 
of individuals (candidate solutions) undergoes evolution following Darwinian 
principles. Each individual is a set of patterns V represented as a vector of 
length equal to the number of shops to serve (nShops), 

-P = (Pl,P2, • • • ,PnShops) 

and whose components pi, 1 < < 2^^ — 1, i G {1 . . .nShops} are integers 
representing a particular delivery pattern, where d is the number of days in the 
working week. 

The relationship between patterns and delivery days is made at bit level. 
Each pattern is formed by d bits and each day corresponds to one bit: 1 means 
that the store is visited that day and that it is not. In our case the working 
week has 5 days (Monday to Friday) so d = 5. Hence patterns are coded by 
the rightmost 5 bits of the integer value. For instance pattern 21, i.e. 10101 in 
binary, corresponds to deliveries on Monday, Wednesday and Friday. 

The total number of possible patterns is 2"^ — 1; in our case 2^ — 1 = 31, so 
patterns range from 1 to 31 or, alternatively, from 00001 to 11111 (obviously 
pattern 00000, i.e. not delivering any day, is not admissible). However, as was 
stated in Subsection 12.21 not all patterns are suitable for all shops. Hence, pi 
must be contained in the set of admissible patterns, pi e Vadm- The elements 
in Vadm are given in Table [21 

Table 2: The admissible patterns and their respective frequencies. The 1 repre- 
sents that the shop is served on that day, the that it is not. As a consequence 
of the business logic, we will only consider 11 patterns out of the 31 that are 
possible. 



Pattern Id. 


Frequency (days) 


Mon 


Tues 


Wed 


Thurs 


Fri 


5 


2 








1 





1 


9 


2 





1 








1 


10 


2 





1 





1 





11 


3 





1 





1 


1 


13 


3 





1 


1 





1 


17 


2 


1 











1 


18 


2 


1 








1 





21 


3 


1 





1 





1 


23 


4 


1 





1 


1 


1 


29 


4 


1 


1 


1 





1 


31 


5 


1 


1 


1 


1 


1 



A point to note is that although the binary representation of the patterns is 
convenient in order to figure out what delivery days are associated to a pattern 
and also for calculating its corresponding frequency, in the genetic algorithm we 
will be considering the integer values only. 
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Algorithm 1 Evaluation function. 



Procedure Evaluate 
input: 

Chromosome {pi, . . . ,p„shops}. 
problem data tables 
costPerKm 
output: Fitness 

[Calculate Inventory cost] 
InventoryCost = 
for i = 1 to nShops 

Look up frequency fi for pattern pi 
Look up cost Ci for shop i and frequency fi 
InventoryCost+ = 
[Calculate transportation cost] 
day = Monday; 
totalDistance = 
repeat 

Identify shops to be served on day 
Run VRP solver to get dayDistance 
total Distance+ = dayDistance 
day + + 
until day = lastW orkingDay; 
TransportCost = totalDistance * costPerKm 
if multiobjective 

return InventoryCost, TransportCost 
[Calculate total cost] 

TotalCost = InventoryCost + TransportCost 
return TotalCost 
end procedure; 



The details of the top level evolutionary algorithm are given in Table [3l The 
pseudo-code for the evaluation function is given in Algorithm [1] 



4 Lower level: Solving the VRP 

The calculation of the fitness of an individual requires an algorithm to solve the 
VRP (a VRPsolver). In this work three algorithms have been tested for this 
purpose, namely: 

CWLS, the classical Clarke and Wright's algorithm ( Clarke and Wrightl , 



19641 ) enhanced with local search, which is presented in subsection 14. II 



ACQ , a bioinspired ant colony optimisation algorithm ( Dorigo and Stutzl3 . 



20041 ) ■ described in subsection 14.21 and 
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Table 3: Configuration of the top-level evolutionary algorithm employed (both 
mono and multiobjective). 



Encoding 


The gene i represents the pattern for shop i. 




The chromosome length is equal to the number of shops 




(nShops). 


Selection 


Tournament in 2 steps. To select each parent, we take 




tSize individuals chosen randomly and select the best. 




For the single objective algorithm the best 10 individuals 




of each generation are preserved as the elite. 


Evolutionary operators 


2 point crossover and 1-point mutation. 




The mutation operator changes the pattern for 1 shop 




in the chromosome. 


Termination criterion 


Terminate when the total number of generations (in- 




cluding the initial one) equals 100. 


Fixed parameters 


Population size, popSize = 100 




Tournament size, tSize = 2 




Mutation probability, pM ^ 0.2 




Crossover probability, pC = 1 



CWTS, tabu search (jGlover and Kochenbergerj . 120021 ) seeded with a solu 



tion obtained by Clarke and Wright algorithm, as described in subsection 



In ( Esparcia- Alcazar et al. . 2009t ) we also employed a fourth algorithm as 



VRPsolver, evolutionary computation. However, the results obtained there were 
not very promising, which is why it has been omitted in this study. 

4.1 Clarke and Wright's algorithm 

Clarke and Wright's algorithm ([Clarke and Wrig htl. Il964h is based on the con- 



cept of saving^ which is the reduction in the traveled length achieved when com- 
bining two routes. We employed the parallel version of the algorithm, which 
works with all routes simultaneously. 

Due to the fact that the solutions generated by the C&W algorithm are not 
guaranteed to be locally optimal with respect to simple neighbourhood defini- 
tions, it is almost always profitable to apply local search to attempt and improve 
each constructed solution. For this purpose we designed a simple and fast local 



9 



search method, which consists on performing 2-interchanges on the solution ob- 
tained by the C&W algorithm. Every possible pair of shops is exchanged, first 
between shops in the same route and then between shops in different routes. If 
at any time an invalid route is generated (because the restrictions on time or 
capacity are violated) the depot is inserted where required in the route. The 
best neighbour solution will be the one with a lower associated transport cost. 

This combination of C&W's algorithm with local search is what we have 
termed CWLS, the pseudo-code for both can be found in Algorithms S] and El 



4.2 Ant Colony Optimisation 



Some ant systems have been applied to the VRP (see for instance (jColtorti and Rizzoli 



2007t iGendreau et al.l . |2003 )) with various degrees of success. Ant algorithms 



are derived from the obs ervation of the self-organized behavior of real ants 
( Dorigo and Stutzld . I200I . The main idea is that artificial agents can imitate 



this behavior and collaborate to solve computational problems by using different 
aspects of ants' behavior. One of the most successful examples of ant algorithms 
is known as "ant colony optimisation", or AGO, which is based on the use of 
pheromones, chemical products dropped by ants when they are moving. 

Each artificial ant builds a solution by choosing probabilistically the next 
node to move to among those it has not visited yet. The choice is biased by 
pheromone trails previously deposited on the graph by other ants and some 
heuristic function. Also, each ant is given a limited form of memory in which 
it can store the partial path it has followed so far, as well as the cost of the 
links it has traversed. This, toge ther with determini stic backward moves, helps 



avoiding the formation of loops (jDorigo and Stutzld . 12004). 



In this work we employed a variant of AGO described by IXiang-pei et al 
(12006 ") which differs from the original AGO algorithm in three aspects: (1) the 
way the pheromone matrix is updated, (2) the transition function and (3) that 
A-interchanges are used instead of local search. 

Two ways of updating the pheromone matrix r are defined: local updating 
and a posteriori updating (i.e. taking place after all ants have built their 
solutions). The former consists of adding l/dij to each element r^.j of the 
pheromone matrix, with dij being the distance between shops i and j. The 
latter is given by Equation ID Here the best path built in iteration t receives a 
reinforcement while the worst path is reset to the initial pheromone value, tq. 



/ , i-p \ if (i,j) is an arc visited by the best ant 

^i.jyP ~^ minCostt ' :„:i„.„4-: -i 



To 



in iteration t 

if (i,j) is an arc visited by the worst ant 
in iteration t 



(4) 

where minCostt is the minimum cost obtained in iteration t and the value 
of p is automatically corrected in each iteration, as follows. 
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^ r 0.95pf_i if 0.95pt_i < pmin 
^* \ Pynin Otherwise ^ 

The second difference with respect to the original AGO is the transition 
function. In our AGO an ant located at shop i will select as its next shop j the 
one given by Equation [6l with probability pt , 

j = argmaxs[Tij]°'[m^jf[p,i.j]'' (6) 

where S is the list of shops not visited yet, r is the pheromone matrix, 77 is 
the heuristic function, 

1 

and finally, 

— di,o + rfoj ~ dij 

corresponds to the concept of saving used in Glarke & Wright's algorithm, 
with shop being the depoo Alternatively, the next shop j will be uniformly 
selected at random from S. 

The parameters a, (3 and 7 (whose sum does not necessarily equal 1) mea- 
sure the relative importance of each component. The probability value pt is 
dynamically adjusted at runtime in a similar way as for pt, following Equation 

m 

^ r 0.95pt_i if 0.95 pt-l > Pmin /-X 

^* \ Pmin Otherwise 

Finally, instead of using local search in the closest neighbours as in con- 
ventional AGO, we use A-interchanges, a concept borrowed from the GWLS 
algorithm described earlier. 

The pseudo-code for the AGO algorithm employed here is given in Algorithm 
[6l its transition function is shown in Algorithm [7] and the parameters used are 
given in Table ID 



4.3 Tabu Search 

Tabu Search (TS) is a mctaheuristic introduced by Glover and Kochen berger in 

order to allow Local Search (LS) methods to overcome local optima (Gl over and Kochenberger . 

The basic principle of TS is to pursue LS whenever a local optimum is 
found by allowing non-improving moves; cycling back to previously visited so- 
lutions is prevented by the use of memories, called tabu lists 
which last for a period given by their tabu tenure. The main TS loop is given 
in Algorithm [51 

°The transition function used in l|Xiang-Dei et al.L 12006^ also takes into account the time 
window of each customer; wc have skipped this because we do not consider time windows in 
our problem. 



GendreaiJ . [l999h 
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Table 4: Par ameters for ACQ . The values were obtained by trial and error, as 
explained in ( Martmez-Garcia , 2008h . Notice that a, (3 and 7 do not necessarily 
add up to 1. 



Number of Iterations 
Number of Ants 
Transition function 



Pheromones 



niters = 50 
nAnts = 25 
Probability, pt 

Initial value, po — 0.8 

Minimum value, Pmin = 0.1 
Weights: 

Pheromone: a — 0.2 

Heuristics: (3 — 0.8 

Savings: 7 — 0.3 
Initial value, tq = 0.5 
Update factor, p 

Initial value, po — 1 

Minimum value, p„iin = 0.1 



To obtain the best neighbour of the current solution we must move in the 
solution's neighbourhood, avoiding moving into older solutions and returning 
the best of all new solutions. This solution may be worse than the current 
solution. For each solution we must generate all possible and valid neighbours 
whose generating moves are not tabu. If a new best neighbour is created, the 
movement is inserted into the Tabu List with the maximum tenure. This move- 
ment is kept in the list until its tenure is over. The tenure can be a fixed or 
variable number of iterations. 

The way to obtain the best neighbour is described in Algorithm HI Table [5] 
lists the relevant information about the algorithm. 



In order to improve over the TS algorithm used in (jEsparcia- Alcazar et al 



200i), we considered seeding the TS with a good solution. The solution chosen 



as a starting point was one obtained by the C&W algorithm. Statistical analysis 
carried ouiLO shows that the seeded algorithm obtains significantly better results 
than when the initial solution is obtained at random. For this reason, in the 
rest of this work we have employed this improved version, which we have termed 
CWTS. 



^Mann- Whitney test on the relative percentage deviation (RPD), given by Eqn of the 
best results of 10 runs per algorithm over the 8 instances of groups A and B, see table El 
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Table 5: Configuration for CWTS algoritlrm. The parameters were tuned heuris- 
tically, after several test launches. The chosen values (small value for the Tabu 
Tenure and large number of iterations without solution improvement) sacrifice 
execution time in order to obtain better solutions. 



Termination criterion 20 iterations without improvement 



5 Experiments 

This section is devoted to present the data we have used in the problem (sub- 
section [231 and the experimental procedure we have followed (in the next sub- 
section, 15. 2p . 

5.1 Problem data 

As explained above, we will employ a number of geographical layouts available 
on the web. We have selected our instances so as to achieve the maximum 
representation on three categories: 

• size, given by the number of shops, 

• distribution. We consider two kinds of distributions: uniform and in 
clusters, corresponding to shops that are scattered more or less uniformly 
on the map or grouped in clusters and, 

• eccentricity. This represents the distance between the depot and the 
geographical centre of the distribution of shops. The coordinates of the 
geographical centre are calculated as follows: 



An instance with low eccentricity (in practise, less than 25) would have 
the depot centered in the middle of the shops while in another with high 
eccentricity (above 40) most shops would be located on one side of the 
depot. 

We chose ten instances with different levels of each category, see Table [HI 



Initial solution 



Tabu tenure 



Possible moves 



A solution obtained by Clarke & Wright's algorithm 
Swap shop i with shop j, both in same route 
Swap shop i with shop j, in different routes 
Create new route with shop i only 
12 iterations 




nShops 



13 



Table 6: Problem instances used in the experiments and 

their characteristics. The first nine have been taken from 

http://branchandcut.org/VRP/data/ and the last one from 
http : //www. f ernuni-hagen .de/WINF/touren/ inhalte/probinst .htm. 



ID 


Instance 


Distribution 


nShops 


Eccentricity 


A32 


A-n32-k5.vrp 


uniform 


31 


47.4 


A33 


A-n33-k5.vrp 


uniform 


32 


20.2 


A69 


A-n69-k9.vrp 


uniform 


68 


15.3 


A80 


A-n80-kl0.vrp 


uniform 


79 


63.4 


B35 


B-n35-k5.vrp 


clusters 


34 


60.5 


B45 


B-n45-k5.vrp 


clusters 


44 


16.6 


B67 


B-n67-kl0.vrp 


clusters 


66 


19.9 


B68 


B-n68-k9.vrp 


clusters 


67 


49.2 


PlOO 


P-nl01-k4.vrp 


uniform 


100 


1.59 


X200 


cl_2_l.txt 


clusters 


200 


8.15 



It must be noted that we are only using the spatial location and not other 
restrictions given in the bibliography, such as the number of vehicles or the shop 
demand values. As pointed out earlier, a main characteristic of our problem is 
that the latter is a function of the delivery frequency, so we had to use our own 
values for the demands. 

We also added a list of admissible patterns, which are given in Table [21 and 
inventory costs, an example of which is given in Table [TJ The inventory, demand 
and admissible patterns data were obtained from Druni SA, a major regional 
Spanish drugstore chain. 

Finally, we have used the vehicle data given in Table [71 



Table 7: Data for Vehicle Routing Problem 



Vehicle capacity 
Transportation cost 
Average speed 
Unloading time 
MELximum working time 



12 roll containers 

0.6 €/Km 

60 km/h 

15 min 

8h 
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5.2 Experimental procedure 



Our aim is, on the one hand, to evaluate whether the multiobjective approach 
yields better results than the single objective one emplo yed in previous work. On 



the o ther, we aim to verify if the conclusions reached in (jEsparcia- Alcazar et al 



I2OO9I) with regard to the best algorithm to use within EVITA still hold after 
the improvements carried out in the TS and AGO algorithms. 

For this purpose, we have tested the selected VRP algorithms described 
above on each one of the ten instances selected, each with a different geographic 
layout. We performed 10 runs per VRP algorithm and instance with a termi- 
nation criterion in all cases of 100 generations. The motivation for such a small 
number of runs is the high computational expense of some of the instance- 
algorithm combinations: the running times ranged from several minutes to sev- 
eral dayf[f| depending on the algorithm and the size of the instance. 

The results were evaluated on two fronts: quantitatively for the total costs 
obtained, and qualitatively for the computational time taken in the runs. The 
latter is important when considering a possible commercial application of the 
EVITA methodology. 



6 Results and analysis 

We carried out Kruskal-Wallis tests for the total costs yielded by the best 
individuals for all runs and problem instances, both in single and multiobjective. 
In the multiobjective case, we define the best individual as that member of 
the final Pareto front yielding the lowest total cost (as defined for the single 
objective problem) The Kruskal-Wallis test is a non-parametric test for multiple 
comparisons which is suited to the case at hand, in which the number of runs 
performed for each combination of instance and VRP solver is small. This test 
does not require normality or homoskedasticity, which are not guaranteed in 
our case. Figures [U [2] and [3] show the resulting boxplots for the ten instances 
studied. 

The boxplot consists of a box and whisker plot for each algorithm. The box 
has lines at the lower quartile, median, and upper quartile values. The whiskers 
are lines extending from each end of the box to show the extent of the rest of 
the data. Outliers are data with values beyond one standard deviation. 

The conclusions that can be reached after analysis of the tests are as follows: 

• The single objective approach yields the best performance in all instances. 

• Considering separately the multi and single objective runs, in general there 
are no significant differences between OWLS and CWTS, although at first 
sight it would seem that CWLS performs better than CWTS on instances 
A32, A69, A80, PlOO and X200 and vice versa on A33, B35, B45 and B67. 
On instance B68 there are no differences at first sight. 

*The computers employed were PCs with Intel Celeron processor, between 1 and 3GHz, 
between 256 and 512 MB RAM. 
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AC0(1 obj) ACO(2obi) CWLS(1 obj) CWLS(2obi) CWTS(1 obj) CWTS(2obi) 
Optimisation methods 



AC0(1 obj) AC0(2 obj) CWLS(1 obj) CWLS(2 obj) CWTS(1 obj) CWTS(2 obj) 
Optimisation methods 



AC0{1 obj) AC0(2 obj) CWLS(1 obj) CWLS(2 obj) CWTS(1 obj) CWTS(2 obj) 
Optimisation methods 
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AC0(1 Obj) ACO{2obi} CWLS(1 obj) CWLS(2obj) CWTS(1 obj) CWTS(2obj) 
Optimisation methods 



Figure 1: Boxplots of total costs for instances A32, A33, B45 and B68. 
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Optimisation methods 
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AC0{1 obj) AC0(2 obj) CWLS(1 obj) CWLS(2 obj) CWTS(1 obj) CWTS(2 obj) 
Optimisation methods 



Figure 2: Boxplots of total costs for instances A69, A80, PlOO and X200. 



• AGO is significantly worse in all cases except B35. Furthermore, it is the 
method that scales worse, getting worse results as the size of the problem 
increases. 

• The case of instance B35 is unique in the sense that monoobjective CWTS 
is significantly better than CWLS and does not differ significantly from 
AGO, both for single and multiobjective. 



AC0(1 obj) AC0(2 obj) CWLS(1 obi) CWLS(2 obi) CWTS(1 obj) CWTS(2 obj) 
Optimisation mettiods 




AC0(1 obi) AC0(2 obj) CWLS(1 obi) CWLS(2 obi) CWTS(1 obi) CWTS(2 obj) 
Optimisation metiiods 



Figure 3: Boxplots of total costs for instances B35 and B67 



In order to be able to compare results between the different instances we 
normalised the fitness values by defining the relative percentage deviation, RPD, 
given by the following expression: 

RPD = - Z^^^^^^"-^" X 100 (8) 
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Figure 4: Boxplots of RPD of total costs for groups clusters (top) and uniform 
(bottom) 



where fitness is the fitness value obtained by an algorithm configuration on 
a given instance. The RPD is, therefore, the average percentage increase over 
the lower bound for each instance, fitnessmin- In our case, the lower bound is 
the best result obtained for that instance across all algorithm configurations. 

With the RPD results of all the runs for all VHP solvers we ran the tests 
again; the results are shown in Figure [4] split into two groups: uniform distri- 
bution of shops and distribution in clusters. The conclusions in this case are 
similar. In both groups CWLS and CWTS perform better than AGO and, at 
first sight, CWLS is better than CWTS for group uniform and vice versa for 
group clusters. Further, considering each VHP solver separately, the single 
objective approach is better than the multiobjective one. 

Figure [5] portrays the comparison between methods when the costs of trans- 
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All Instances 



AC0(1 obi) AC0(2 obi) CWLS(1 obj) CWLS(2 obj) CWTS(1 obi) CWTS(2 obi) 
Optimisation methods 



AC0(1 obi) AC0(2 obj) CWLS(1 obj) CWLS(2 obj) CWTS(1 obj) CWTS(2 obj) 
Optimisation methods 



Figure 5: Boxplots for the transport (top) and inventory costs (bottom) 
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port and inventory are considered separately. Here it can be observed that there 
are no significant difi^erences between methods if only inv entory costs are taken 



into account. This differs from the conclusions obtained in flEsparcia- Alcazar et al 



2006a(), where the choice of VRPsolver influenced the inventory cost results. In 



that work, however, the algorithms employed were suboptimal compared to the 
ones used here. So, it could be concluded that the choice of VRPsolver does not 
have an influence on the inventory costs provided a "good enough" algorithm is 
chosen. 

The big difference lies in the transport costs, which is where AGO clearly 
shows its inferiority, especially in the multiobjective approach. 

Regarding the computational time, the results clearly favour GWLS over all 
other VRP solvers. In general, when employing GWLS the time for a whole 
run took approximately the same as that of a single generation in when using 
GWTS or AGO. This is a point in favour of GWLS when considering a potential 
commercial application. 



6.1 Pareto fronts 

In this section we study the Pareto fronts obtained in the multiobjective ap- 
proach in two instances of the problem, namely B35 and A32. The former has 
been chosen because of its uncharacteristic behaviour (as we have seen, in this 
instance AGO performs better than the other VRP solvers), and the latter be- 
cause it is of a similar size. Pareto fronts for both instances are shown in Figure 
[6l In Table [8] we show the number of non-dominated solutions vs. the number 
of different individuals for each algorithm. 



Table 8: Final population characteristics in multi-objective approach for one 
run of instances A32 and B35 using GWLS, AGO and GWTS as VRP-solvers. 
Population size is 100 individuals. 



Instance A32 




No. of different individuals 


Pareto front size 


Ratio 


CWLS 

AGO 

CWTS 


7 

100 
24 


7 
19 

24 


1 

0.19 
1 


Instance B35 




No. of different individuals 


Pareto front size 


Ratio 


CWLS 

AGO 

CWTS 


18 
100 
36 


18 
8 

26 


1 

0.08 
0.72 



Comparing the number of non-dominated solutions in the final generation for 
each instance, we observe that the Pareto front yielded by the best performing 
VRP solver (GWLS in A32 and AGO in B35) has always a smaller size, whilst 
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Figure 6: Pareto fronts in multi-objective approach for one run of instances A32 
(left) and B35 (right) using CWLS, ACO and CWTS as VRP-solvers. 



the one corresponding to the worst performing VRP solver has a bigger size 
(ACO in A32 and CWLS in B35). The figures also hint at the possibility that 
using CWLS or ACO as the VRP solver causes the Pareto front to converge to 
a very reduced set of solutions, while when using CWTS more non-dominated 
solutions are preserved. 

On the other hand the ratio of repeated individuals over the population size 
is very high when using CWTS or CWLS, whilst in ACO there are no repeated 
individuals. It is possible that algorithm ACO has a slower rate of convergence 
than the other two in most of instances (except in cases such as B35), so it 
requires more generations to obtain equivalent results. 

At this point it would be of interest to measure the quality of t he Pareto front 
using one of the metrics that have been proposed for this purpose (jCoello Coellol 
20051 ). However, many of them (e.g. the error ratio or the generational distance) 



assume that knowledge exists on the actual Pareto front, which is not the case 
here. Other metrics measure the distribution of solutions on the Pareto front 
by evaluating the variance of neighboring solutions. An example of this is the 
spacing, SP, which measures the relative distances between the members of 
Pareto front; a value of SP = means that members of the Pareto front are 
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equispaced. The spacing is given by the following equation: 



SP=^ 



\ 



1 " 



1=1 

where n is the number of non-dominated solutions found, the distance di is 
given by 

d, = minj{\fl{x) - f({x)\ + If^ix) - fi{x)\), i,j = 1, ...,n 

where f^{-) is the fitness of point k on objective N, x is the generation 
number and d is the mean of all di . 

The values obtained are given in Tabled From these we can see that the 
best values of the metric (i.e. the lowest spacing) are obtained by CWTS; 
however, from previous analysis we know that it is the other two VRP solvers 
that perform better: AGO for B35 and CWLS for A32. We can hence conclude 
that this metric is not very meaningful for the purposes of our problem. 



Table 9: Values of the spacing for the Pareto fronts generated by all VRP solvers 



VRP solver 


Problem instance 


A32 


B35 


CWLS 

AGO 

CWTS 


17.605 
12.823 
8.394 


23.937 
68.027 
17.858 



7 Conclusions and future work 

We have shown how, for the problem presented here, the multiobjective ap- 
proach does not yield any advantage over the single objective one. This could 
be explained by the fact that inventory costs are well above the transport costs. 
Given that the multiobjective approach does not prefer one objective over the 
other, it can happen that there are solutions for which transport costs are very 
low, but still that does not compensate for high inventory costs. 

Further, we have shown how a classical algorithm such as Clarke and Wright's, 
enhanced with local search, can be the best choice in the context of the Inven- 
tory and Transportation Problem, both in terms of the quality of the solutions 
obtained and the computational time necessary to achieve them. The power of 
other algorithms known to perform well in the context of VRP, such as AGO 
and TS, does not grant a good performance for the joint inventory and trans- 
portation problem. In general, using a global optimisation algorithm such as 
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evolutionary computation jointly with a heuristic method adapted to the pro- 
blem at hand, such as CWLS, yields the best results, so this is no surprise. 

It could be argued that both TS and AGO require a finer tuning of their 
parameters in order to give an adequate performance than what was achieved 
here. This, however, could be interpreted as a disadvantage of their application 
to a variety of problem configurations and in a commercial context. 

Special attention should be given to the case of instance B35, since it is the 
only one for which AGO yields better results than the remaining VRP solvers 
(with the exception of single objective TS) and, oddly, this happens in the 
multiobjective case. The geographical layout of this instance is shown in Figure 
[71 the high eccentricity of the distribution can be seen (the depot is on one side 
of the shops) compared to the small number of shops. Instance A80 also has a 
similar value of eccentricity, but for a higher number of shops. Perhaps here lies 
the explanation of why AGO works better in the former but not in the latter; 
clarifying this point is left for future work. In any case, we can conclude that 
although the GWLS with single objective approach performs better in general, it 
is nonetheless interesting to have a tool that can provide the possibility of using 
other VRP solvers and a multiobjective approach in order to handle special 
cases, such as B35. 
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Figure 7: Geographical layout for instance B35 



It must be noted that this result is subject to the specificities of the problem 
data, i.e. the fact that in this case the inventory cost greatly outweighs the 
cost of the transport. As future work we must consider a case in which the 
products moved are cheaper and hence the inventory cost is more in a par with 
the transport cost. 
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Appendix I: NSGA-II 

NSGA-II (|Deb et all l200l l2000l) is an non-elitist multiobjective evolutionary 
algorithm (MOEA) which was developed in order to overcome the problems of 
previous MOEAs, such as the high computational complexity of sorting non 
dominated solutions. These algorithms have two common features: assigning 
fitness to population members based on nondominated sorting and preserving 
the diversity among solutions of the same nondominated front. 

NSGA-II works as follows: Initially, a random parent population Pq is cre- 
ated, with size N. The population is sorted based on nondominance. Each 
solution is assigned a fitness (or rank) equal to its nondomination level (with 1 
being the best level). Initially, the usual binary tournament selection, recombi- 
nation, and mutation operators are used to create an offspring population Qo of 
size N. For the remaining generations t we do the following: First, a population 
Rt of size 2N is formed as the union of Pt and Qt ^ind sorted according to non- 
domination in a number of fronts F. Next, a new population Pt+i of size N is 
created by selecting individuals from Rt in order of nondominance (i.e. ordered 
individuals from Fi are chosen first, then ordered individuals from F2 and so 
on until the number of individuals belonging to Pt+i is N). 

The advantages of NSGA-II with respect to previous MOEAs are the fast 
sorting of nondominated individuals and the preservation of diversity. 

Nondominated sorting. First, for each solution p we calculate two entities: 
domination count, Up, i.e. the number of solutions which dominate p, and a list 
of solutions that p dominates. Bp. All solutions in the first nondominated front 
will have their Up as zero. Now, for each solution with Up ~ 0, we visit each 
member q in its Sp and reduce its domination count by one. In doing so, if for 
any member the domination count becomes zero, we place it in a separate list 
Q. These members belong to the second nondominated front. This procedure is 
continued with each member of Q and the third front is identified. This process 
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continues until all fronts are identified. The code for this operation can be found 
in Algorithm [21 



Algorithm 2 Nondominated sorting function for NSGAII. 



Procedure fastNonDominatedSort() 
input: population P 
output: listOfFronts F 
For each p £ P 

Sp = 
Up = 

For each q a P 

if p dominates q 

Adds q to Sp 
else if q dominates p 

Increments rip 

if rtp = 

p.rank = 1 

Adds p to F\ 
/ = 1 

while Ff 7^ 

Q = % 

for each p Ff 

for each q Sp 

Decrements Uq 
if rtg = 

q.rank = 1 
Adds q to Q 

Ff = Q 
Increments / 

return F 
end; 



Diversity preservation. To get an estimate of the density of solutions sur- 
rounding a particular solution in the population, we calculate the average dis- 
tance of two points on either side of this point along each of the objectives. 
This distance is called the crowding distance, and it is calculated as shown in 
Algorithm [3l Moreover, a crowded- comparison operator,<n, is used in order 
to guide the selection process at the various stages of the algorithm towards a 
uniformly spread-out Pareto-optimal front. 

In the selection process, given two solutions with different nondomination 
ranks the one with the lower (better) rank will be preferred. Otherwise, if both 
solutions belong to the same front, then one located in a less crowded region 
will be preferred. 

For a complete description of NSGA-II, see (jPeb et all . I200ll2000( ). 
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Algorithm 3 Crowding distance assignment for NSGA-II Algorithm. 



Procedure crowdingDistance() 
input: population pop[popSize] 

for i = L.popSize 

pop[i]. distance = 
for each objective m 

Sort pop using m 

pop[l]. distance = pop\popSize]. distance = oo 
for i = 2..popSize — 1 

pop[i].distance = pop[i].distance + {pop[i + l].m — pop[i — l].m)/{fj^°'^ — /^™) 

end; 



Appendix II: Algorithms for solving the VRP 

Here wc describe the algorithms employed as VRPsolvers: C&W's algorithm, 
Local Search, AGO and TS, plus other additional functions. 
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Algorithm 4 Clarke & Wright's algorithm (C&W). 



Algorithm CWLS 

input: shops[nShops], depot 

output: routes 

[Build nShops initial routes with one sHiop only] 

for 1=1... nShops 

Ti = {depot, shops[i], depot) 
[Calculate savings] 

Calculate Sij for each pair shops[i], shopsp] 

^^^shops[i]depot ~^ ^^^depotshopslj] ^^ishops[i]shops[j] 

[Best union ] 
repeat 

Let ria, be the route containing shops[i] 
Let rjt be the route containing shops[j] 
if 

shops[i*] is the last shop in 
and shopsp*] is the first shop in r^* 
and the combination is feasible 
then combine rj* and rj* 
delete Sit^*; 
until there are no more savings to consider; 
return routes 

end; 
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Algorithm 5 Local Search (LS) algorithm. 



Algorithm LocalSearch 
input: initialSolution 
[Initialisation] 

best <— initialSolution 
costBest = costTemp 
[Improving each route] 
for each r in routes(initialSolution) 

for k= 1 . . M ax(sizeOf ( r) , num berOf Neighbors) 
select si, s2 random shops in r 

tempSolution <— lnterchangeShops(initialSolution,sl, s2) 
if (calculateCost(tempSolution) < costBest) 
best -t— tempSolution 
costBest = calculateCost(best) 
[Improving pairs of routes] 
for rl =l..routes(initialSolution) 

for r2=rl+l..routes(initialSolution) 
select si random shop in rl 
select s2 random shop in r2 

solTemp <— lnterchangeShops(initialSolution,sl, s2) 
if (calculateCost(tempSolution) < costBest) 

best <— tempSolution 

costBest = calculateCost(best) 

return best 
end algorithm; 
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Algorithm 6 AGO - Main loop. 



Algorithm AGO 

input: shops[nShops], numberOflterations 
output: routes 
[Initialise values] 

routes <— validRandomSolution 
globalCost <— calculateCost(routes) 
[Main loop] 

for it— 1.. numberOflterations 

[Looking for a solution for each ant] 
for each ant i 

Place ant i at depot 
i.shopsNotVisited shops 
i. solution ^ 

i. solution reset cost, tinne and demand 
reset load 

while there are ants i for which I.shopsNotVisited do 
for each ant i for which I.shopsNotVisited 7^ 
nextShopToVisit <— TransitionFunction() 
Update cost 

if nextShopToVisit == depot then reset time and load 
else 

Update time, demand and load 
Update i.shopsNotVisited, i. solution 
Place ant i at nextShopToVisit 
[Interchanges] 

for each ant i do A-interchanges in i. solution 
[Update pheromone matrix with local solution] 
Find worstPathlnlteration, bestPathlnlteration 
Reinforce bestPathlnlteration in pheromone matrix 
Reset worstPathlnlteration in pheromone matrix 
[Update global solution] 

if globalCost > cost(bestPathlnlteration) 
routes -f— bestPathlnlteration 
globalCost -f— cost(bestPathlnlteration) 
[Update configuration] 

Update algorithm parameters: p, p 
return routes 
end algorithm; 
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Algorithm 7 AGO - Transition function. 



Algorithm Transition(shopsNotVisitedYet:List(shop),currentShop:shop, 
remainingTime:float,remainingLoad:float):nextShop 

j = currentShop 
if (g <pt) 

j ^^3^^-^shopsNotVisitedYet[^currentShop,j] [VcurrentShop,j]^\McurrentShop,j\ 

else 

j = randomFrom{shopsNotVisitedYet) 
time = TimeO f {currentShop, j) + TimeOf{j, Depot) 
if {time < remainingTime AND DemandOf{j) < remaining Load ) 

nextShop = j 

else 

nextShop = Depot 
return nextShop 
end algorithm; 



Algorithm 8 Tabu Search (TS) algorithm - Main loop. 

Algorithm TS 

[Initialisation] 

currentSolution ^ initialSolution 

currentSolutionCost <— calculateCost(currentSolution) 
[Main loop] 

while (iterations < MAX-ITERATIONS) 

bestNeighbour <— getBestNeighbour (currentSolution, tabuList) 
bestNeighbourCost «— calculateCost(bestNeighbour) 
currentSolution <— bestNeighbour 
currentSolutionCost <— bestNeighbourCost 
if (currentSolutionCost. isBetterThan(bestsolutionCost)) 
bestSolution <— currentSolution 
bestSolutionCost currentSolutionCost 
iterations <— 

else 

iterations <— iterations + 1 
tabuList <— updateTenure 
return bestSolution 
end algorithm; 
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Algorithm 9 Taljii St^arcli - B(>sl uciglilxjur Mlgoritlini. 



Algorithm bestNeighbour 
[Initialisation] 

moved ^ false 
moves <— getAIIMoves 
theBestNeighbour <— currentSolution 
theBestNeighbourCost <— oo 
neighbourCost <— oo 
[Main loop] 

for i=l:moves.length 

move <— moves[i] 

neighbour currentSolution 

neighbour <— move.operateOn(neighbour) 

neighbourCost •<— calculateCost(neighbour) 

isTabu ^ isTabu(move) 

[Aspiration criteria] 

if (neighbourCost < bestSolutionCost) 
isTabu <— false 

if (neighbourCost < theBestNeighbourCost AIMD NOT isTabu) 
theBestNeighbour ^ neighbour 
theBestNeighbourCost <— neighbourCost 
bestNeighbourMove <— move 
moved <— false 

end for 
[Update tabu list] 

if moved == true 

tabuList.addMove(bestNeighbourMove) 
return theBestNeighbour 
end algorithm; 
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